Fit VBGE models to back-calculated length-at-age

Author

Max Lindmark, Jan Ohlberger, Anna Gårdmark

Published

June 14, 2023

Load libraries

Code
pkgs <- c("tidyverse", "tidylog", "rTPC", "nls.multstart", "broom", "RColorBrewer", "viridis", "minpack.lm", "patchwork")

## minpack.lm needed if using nlsLM()
if(length(setdiff(pkgs, rownames(installed.packages()))) > 0){

    install.packages(setdiff(pkgs, rownames(installed.packages())), dependencies = T)
  
  }

invisible(lapply(pkgs, library,character.only = T))
#> ── Attaching core tidyverse packages ──────────────────────── tidyverse 2.0.0 ──
#> ✔ dplyr     1.1.2     ✔ readr     2.1.4
#> ✔ forcats   1.0.0     ✔ stringr   1.5.0
#> ✔ ggplot2   3.4.2     ✔ tibble    3.2.1
#> ✔ lubridate 1.9.2     ✔ tidyr     1.3.0
#> ✔ purrr     1.0.1     
#> ── Conflicts ────────────────────────────────────────── tidyverse_conflicts() ──
#> ✖ dplyr::filter() masks stats::filter()
#> ✖ dplyr::lag()    masks stats::lag()
#> ℹ Use the conflicted package (<http://conflicted.r-lib.org/>) to force all conflicts to become errors
#> 
#> Attaching package: 'tidylog'
#> 
#> 
#> The following objects are masked from 'package:dplyr':
#> 
#>     add_count, add_tally, anti_join, count, distinct, distinct_all,
#>     distinct_at, distinct_if, filter, filter_all, filter_at, filter_if,
#>     full_join, group_by, group_by_all, group_by_at, group_by_if,
#>     inner_join, left_join, mutate, mutate_all, mutate_at, mutate_if,
#>     relocate, rename, rename_all, rename_at, rename_if, rename_with,
#>     right_join, sample_frac, sample_n, select, select_all, select_at,
#>     select_if, semi_join, slice, slice_head, slice_max, slice_min,
#>     slice_sample, slice_tail, summarise, summarise_all, summarise_at,
#>     summarise_if, summarize, summarize_all, summarize_at, summarize_if,
#>     tally, top_frac, top_n, transmute, transmute_all, transmute_at,
#>     transmute_if, ungroup
#> 
#> 
#> The following objects are masked from 'package:tidyr':
#> 
#>     drop_na, fill, gather, pivot_longer, pivot_wider, replace_na,
#>     spread, uncount
#> 
#> 
#> The following object is masked from 'package:stats':
#> 
#>     filter
#> 
#> 
#> Loading required package: viridisLite

# devtools::install_github("seananderson/ggsidekick") ## not on CRAN 
library(ggsidekick)
theme_set(theme_sleek())

# Check the temperature model script! This is the order based on mean temperature, which makes sense for the sharpe schoolfield plot, and therefore we might keep it across plots
order <- c("SI_HA", "BT", "TH", "SI_EK", "FM", "VN", "JM", "MU", "FB", "BS", "HO", "RA")

nareas <- length(unique(order))
colors <- colorRampPalette(brewer.pal(name = "RdYlBu", n = 10))(nareas)

# Load functions
home <- here::here()
fxn <- list.files(paste0(home, "/R/functions"))
invisible(sapply(FUN = source, paste0(home, "/R/functions/", fxn)))

Read and trim data

Code
d <- #read_csv(paste0(home, "/data/for_analysis/dat.csv")) |> 
  read_csv("https://raw.githubusercontent.com/maxlindmark/perch-growth/master/data/for_analysis/dat.csv") |> 
  filter(age_ring == "Y") # use only length-at-age by filtering on age_ring
#> Rows: 364546 Columns: 11
#> ── Column specification ────────────────────────────────────────────────────────
#> Delimiter: ","
#> chr (5): age_ring, area, gear, ID, sex
#> dbl (6): length_mm, age_bc, age_catch, catch_year, cohort, final_length
#> 
#> ℹ Use `spec()` to retrieve the full column specification for this data.
#> ℹ Specify the column types or set `show_col_types = FALSE` to quiet this message.
#> filter: removed 28,878 rows (8%), 335,668 rows remaining

# Sample size by area and cohort
ns <- d |> 
  group_by(cohort, area) |> 
  summarise(n = n())
#> group_by: 2 grouping variables (cohort, area)
#> summarise: now 505 rows and 3 columns, one group variable remaining (cohort)

# Minimum number of observations per area and cohort
d$area_cohort <- as.factor(paste(d$area, d$cohort))

d <- d |>
  group_by(area_cohort) |> 
  filter(n() > 100)
#> group_by: one grouping variable (area_cohort)
#> filter (grouped): removed 2,844 rows (1%), 332,824 rows remaining

# Minimum number of observations per area, cohort, and age
d$area_cohort_age <- as.factor(paste(d$area, d$cohort, d$age_bc))

d <- d |>
  group_by(area_cohort_age) |> 
  filter(n() > 10)
#> group_by: one grouping variable (area_cohort_age)
#> filter (grouped): removed 3,974 rows (1%), 328,850 rows remaining

# Minimum number of cohorts in a given area
cnt <- d |>
  group_by(area) |>
  summarise(n = n_distinct(cohort)) |>
  filter(n >= 10)
#> group_by: one grouping variable (area)
#> summarise: now 12 rows and 2 columns, ungrouped
#> filter: no rows removed

d <- d[d$area %in% cnt$area, ]

# Plot cleaned data
ggplot(d, aes(age_bc, length_mm)) +
  geom_jitter(size = 0.1, height = 0, alpha = 0.1) +
  scale_x_continuous(breaks = seq(20)) +
  theme(axis.text.x = element_text(angle = 0)) +
  theme(axis.text = element_text(size = 12), axis.title = element_text(size = 15)) +
  labs(x = "Age", y = "Length (mm)") +
  guides(color = "none") + 
  facet_wrap(~area, scale = "free_x")

Code

# Longitude and latitude attributes for each area
area <- c("BS", "BT", "FB", "FM", "HO", "JM", "MU", "RA", "SI_EK", "SI_HA", "TH", "VN")
nareas <- length(area)
lat <- c(60, 60.4, 60.3, 60.5, 63.7, 58, 59, 65.9, 57.3, 57.4, 56.1, 57.5)
lon <- c(21.5, 18.1, 19.5, 18, 20.9, 16.8, 18.1, 22.3, 16.6, 16.7, 15.9, 16.9)
area_attr <- data.frame(cbind(area = area, lat = lat, lon = lon)) |>
  mutate_at(c("lat","lon"), as.numeric)
#> mutate_at: converted 'lat' from character to double (0 new NA)
#>            converted 'lon' from character to double (0 new NA)

Fit VBGE models

Code
# Get individual growth parameters (functions: VBGF/Gompertz,nls_out,fit_nls)
IVBG <- d |> 
  group_by(ID) |> 
  summarize(k = nls_out(fit_nls(length_mm, age_bc, min_nage = 4, model = "VBGF"))$k,
            k_se = nls_out(fit_nls(length_mm, age_bc, min_nage = 4, model = "VBGF"))$k_se,
            linf = nls_out(fit_nls(length_mm, age_bc, min_nage = 4, model = "VBGF"))$linf,
            linf_se = nls_out(fit_nls(length_mm, age_bc, min_nage = 4, model = "VBGF"))$linf_se)
#> group_by: one grouping variable (ID)
#> summarize: now 91,422 rows and 5 columns, ungrouped

Inspect predictions

Code
IVBG <- IVBG |> drop_na(k) # The NAs are because the model didn't converge or because they were below the threshold age
#> drop_na: removed 50,514 rows (55%), 40,908 rows remaining

IVBG <- IVBG |>
  mutate(k_lwr = k - 1.96*k_se,
         k_upr = k + 1.96*k_se,
         linf_lwr = linf - 1.96*linf_se,
         linf_upr = linf + 1.96*linf_se,
         row_id = row_number())
#> mutate: new variable 'k_lwr' (double) with 40,864 unique values and 0% NA
#>         new variable 'k_upr' (double) with 40,864 unique values and 0% NA
#>         new variable 'linf_lwr' (double) with 40,864 unique values and 0% NA
#>         new variable 'linf_upr' (double) with 40,864 unique values and 0% NA
#>         new variable 'row_id' (integer) with 40,908 unique values and 0% NA

# Plot all K's
IVBG |>
  #filter(row_id < 5000) |>
  ggplot(aes(row_id, k, ymin = k_lwr, ymax = k_upr)) +
  geom_point(alpha = 0.5) +
  geom_errorbar(alpha = 0.5) +
  NULL

Code

# Plot all L_inf's
IVBG |>
  #filter(row_id < 5000) |>
  ggplot(aes(row_id, linf, ymin = linf_lwr, ymax = linf_upr)) +
  geom_point(alpha = 0.5) +
  geom_errorbar(alpha = 0.5) +
  NULL

Code

# We can also consider removing the upper 95% quantile of standard errors (not quantile of K directly)
IVBG |>
  tidylog::mutate(keep = ifelse(k > quantile(k_se, probs = 0.95), "N", "Y")) |>
  #filter(row_id < 10000) |>
  ggplot(aes(row_id, k, ymin = k_lwr, ymax = k_upr, color = keep)) +
  geom_point(alpha = 0.5) +
  facet_wrap(~keep, ncol = 1) +
  geom_errorbar(alpha = 0.5) +
  NULL
#> mutate: new variable 'keep' (character) with 2 unique values and 0% NA

Code

# Add back cohort and area variables
IVBG <- IVBG |> 
  left_join(d |> select(ID, area_cohort)) |> 
  separate(area_cohort, into = c("area", "cohort"), remove = TRUE, sep = " ") |> 
  mutate(cohort = as.numeric(cohort))
#> Adding missing grouping variables: `area_cohort_age`
#> select: dropped 10 variables (length_mm, age_bc, age_catch, age_ring, area, …)
#> Joining with `by = join_by(ID)`
#> left_join: added 2 columns (area_cohort_age, area_cohort)
#> > rows only in x 0
#> > rows only in y (115,512)
#> > matched rows 213,338 (includes duplicates)
#> > =========
#> > rows total 213,338
#> mutate: converted 'cohort' from character to double (0 new NA)

# Summarise and save for sample size
sample_size <- IVBG |>
  group_by(area) |> 
  summarise(n_cohort = length(unique(cohort)),
            min_cohort = min(cohort),
            max_cohort = min(cohort),
            n_individuals = length(unique(ID)),
            n_data_points = n())
#> group_by: one grouping variable (area)
#> summarise: now 12 rows and 6 columns, ungrouped

write_csv(sample_size, paste0(home, "/output/sample_size.csv"))

# Compare how the means and quantiles differ depending on this filtering
IVBG_filter <- IVBG |>
  tidylog::drop_na(k_se) |>
  tidylog::filter(k_se < quantile(k_se, probs = 0.95))
#> drop_na: no rows removed
#> filter: removed 10,670 rows (5%), 202,668 rows remaining

# Summarize growth coefficients by cohort and area
VBG <- IVBG |>
  group_by(cohort, area) |>
  summarize(k_mean = mean(k, na.rm = T),
            k_median = quantile(k, prob = 0.5, na.rm = T),
            linf_median = quantile(linf, prob = 0.5, na.rm = T),
            k_lwr = quantile(k, prob = 0.05, na.rm = T),
            k_upr = quantile(k, prob = 0.95, na.rm = T)) |> 
  ungroup()
#> group_by: 2 grouping variables (cohort, area)
#> summarize: now 382 rows and 7 columns, one group variable remaining (cohort)
#> ungroup: no grouping variables

VBG_filter <- IVBG_filter |>
  group_by(cohort, area) |>
  summarize(k_mean = mean(k, na.rm = T),
            k_median = quantile(k, prob = 0.5, na.rm = T),
            k_lwr = quantile(k, prob = 0.05, na.rm = T),
            k_upr = quantile(k, prob = 0.95, na.rm = T)) |> 
  ungroup()
#> group_by: 2 grouping variables (cohort, area)
#> summarize: now 382 rows and 6 columns, one group variable remaining (cohort)
#> ungroup: no grouping variables

ggplot() +
  geom_ribbon(data = VBG, aes(cohort, k_median, ymin = k_lwr, ymax = k_upr, fill = "All k's"), alpha = 0.4) +
  geom_ribbon(data = VBG_filter, aes(cohort, k_median, ymin = k_lwr, ymax = k_upr, fill = "Filtered"), alpha = 0.4) +
  geom_line(data = VBG, aes(cohort, k_median, color = "All k's")) + 
  geom_line(data = VBG_filter, aes(cohort, k_median, color = "Filtered")) + 
  guides(fill = FALSE) +
  facet_wrap(~area)
#> Warning: The `<scale>` argument of `guides()` cannot be `FALSE`. Use "none" instead as
#> of ggplot2 3.3.4.

Code

ggplot() +
  geom_line(data = VBG, aes(cohort, k_median, color = "All k's")) + 
  geom_line(data = VBG_filter, aes(cohort, k_median, color = "Filtered")) + 
  guides(fill = FALSE) +
  facet_wrap(~area)

Code

# No difference at all really. We should probably just discuss that with this model, achieving biologically reasonable parameter values and a good fit to data are sometimes two different things. In our case, we just want a representative value of the growth (as in time to reach average maximum size in the population) that accounts for the entire growth trajectory of an individual, for each area and cohort.

Add GAM-predicted temperature to growth models

Code
pred_temp <- read_csv(paste0(home, "/output/gam_predicted_temps.csv")) |> 
  rename(cohort = year)
#> Rows: 566 Columns: 4
#> ── Column specification ────────────────────────────────────────────────────────
#> Delimiter: ","
#> chr (2): area, type
#> dbl (2): year, temp
#> 
#> ℹ Use `spec()` to retrieve the full column specification for this data.
#> ℹ Specify the column types or set `show_col_types = FALSE` to quiet this message.
#> rename: renamed one variable (cohort)

VBG <- VBG |> left_join(pred_temp, by = c("area", "cohort"))
#> left_join: added 2 columns (temp, type)
#>            > rows only in x     0
#>            > rows only in y  (184)
#>            > matched rows     382
#>            >                 =====
#>            > rows total       382

# save data for map-plot
cohort_sample_size <- IVBG |>
  group_by(area, cohort) |> 
  summarise(n = n()) # individuals, not samples!
#> group_by: 2 grouping variables (area, cohort)
#> summarise: now 382 rows and 3 columns, one group variable remaining (area)
  
VBG <- left_join(VBG, cohort_sample_size, by = c("cohort", "area"))
#> left_join: added one column (n)
#>            > rows only in x     0
#>            > rows only in y  (  0)
#>            > matched rows     382
#>            >                 =====
#>            > rows total       382

write_csv(VBG, paste0(home, "/output/vbg.csv"))

Plot VBGE fits

Code
ids <- d |> group_by(area) |> 
  sample_n(50)
#> group_by: one grouping variable (area)
#> sample_n (grouped): removed 328,250 rows (>99%), 600 rows remaining

dat <- IVBG %>%
  filter(ID %in% ids$ID)
#> filter: removed 211,573 rows (99%), 1,765 rows remaining

d |>
  ungroup() |> 
  filter(ID %in% ids$ID) |>
  left_join(dat) |> 
  tidylog::mutate(length_mm_pred = linf*(1-exp(-k*age_bc))) |>
  ggplot(aes(age_bc, length_mm, group = ID, color = ID)) +
  geom_jitter(width = 0.2, height = 0, alpha = 0.5, size = 0.6) +
  geom_line(aes(age_bc, length_mm_pred, group = ID, color = ID), inherit.aes = FALSE, alpha = 0.5) + 
  guides(color = "none") +
  scale_color_viridis(discrete = TRUE, name = "Area", option = "cividis") +
  labs(x = "Age (yrs)", y = "Length (mm)") +
  facet_wrap(~factor(area, order))
#> ungroup: no grouping variables
#> filter: removed 326,388 rows (99%), 2,462 rows remaining
#> Joining with `by = join_by(area, cohort, ID, area_cohort_age)`
#> left_join: added 9 columns (k, k_se, linf, linf_se, k_lwr, …)
#> > rows only in x 697
#> > rows only in y ( 0)
#> > matched rows 1,765
#> > =======
#> > rows total 2,462
#> mutate: new variable 'length_mm_pred' (double) with 1,766 unique values and 28%
#> NA
#> Warning: Removed 697 rows containing missing values (`geom_line()`).

Code

ggsave(paste0(home, "/figures/vb_fits.pdf" ), width = 17, height = 17, units = "cm")
#> Warning: Removed 697 rows containing missing values (`geom_line()`).

k <- IVBG |> 
  ggplot(aes(k, color = factor(area, order))) + 
  geom_density(alpha = 0.4, fill = NA) +  
  scale_color_viridis(discrete = TRUE, name = "Area", direction = -1) +
  coord_cartesian(expand = 0) + 
  guides(color = "none") +
  theme(legend.position = c(0.8, 0.8))

linf <- IVBG |> 
  ggplot(aes(linf, color = factor(area, order))) + 
  geom_density(alpha = 0.4, fill = NA) +  
  scale_color_viridis(discrete = TRUE, name = "Area", direction = -1) +
  coord_cartesian(expand = 0, xlim = c(0, 2000)) + 
  theme(legend.position = c(0.9, 0.5))

k / linf

Code

ggsave(paste0(home, "/figures/vb_pars.pdf" ), width = 17, height = 17, units = "cm")

Fit Sharpe-Schoolfield model to K

By area

Code
model <- 'sharpeschoolhigh_1981'

# Get starting values on full dataset for Sharpe-Schoolfield model
dat <- VBG |>
  select(k_median, temp) |>
  rename(rate = k_median)
#> select: dropped 8 variables (cohort, area, k_mean, linf_median, k_lwr, …)
#> rename: renamed one variable (rate)

lower <- get_lower_lims(dat$temp, dat$rate, model_name = model)
upper <- get_upper_lims(dat$temp, dat$rate, model_name = model)
start <- get_start_vals(dat$temp, dat$rate, model_name = model)
  
# Sharpe-Schoolfield model fit to data for each area
preds <- NULL

for(a in unique(VBG$area)) {
  
  # Get data
  dat <- VBG |> 
    filter(area == a) |> 
    select(k_median, temp, area) |> 
    rename(rate = k_median)
  
  # Fit model
  fit <- nls_multstart(
    rate ~ sharpeschoolhigh_1981(temp = temp, r_tref, e, eh, th, tref = 8),
    data = dat,
    iter = c(3, 3, 3, 3),
    start_lower = start*0.5,
    start_upper = start*2,
    lower = lower,
    upper = upper,
    supp_errors = 'Y'
    )
  
  # Make predictions on new data
  new_data <- data.frame(temp = seq(min(dat$temp), max(dat$temp), length.out = 100))
  
  pred <- augment(fit, newdata = new_data) |>
    mutate(area = a)
  
  # Add to general data frame
  preds <- data.frame(rbind(preds, pred))
  
}
#> filter: removed 319 rows (84%), 63 rows remaining
#> select: dropped 7 variables (cohort, k_mean, linf_median, k_lwr, k_upr, …)
#> rename: renamed one variable (rate)
#> mutate: new variable 'area' (character) with one unique value and 0% NA
#> filter: removed 328 rows (86%), 54 rows remaining
#> select: dropped 7 variables (cohort, k_mean, linf_median, k_lwr, k_upr, …)
#> rename: renamed one variable (rate)
#> mutate: new variable 'area' (character) with one unique value and 0% NA
#> filter: removed 341 rows (89%), 41 rows remaining
#> select: dropped 7 variables (cohort, k_mean, linf_median, k_lwr, k_upr, …)
#> rename: renamed one variable (rate)
#> mutate: new variable 'area' (character) with one unique value and 0% NA
#> filter: removed 347 rows (91%), 35 rows remaining
#> select: dropped 7 variables (cohort, k_mean, linf_median, k_lwr, k_upr, …)
#> rename: renamed one variable (rate)
#> mutate: new variable 'area' (character) with one unique value and 0% NA
#> filter: removed 334 rows (87%), 48 rows remaining
#> select: dropped 7 variables (cohort, k_mean, linf_median, k_lwr, k_upr, …)
#> rename: renamed one variable (rate)
#> mutate: new variable 'area' (character) with one unique value and 0% NA
#> filter: removed 344 rows (90%), 38 rows remaining
#> select: dropped 7 variables (cohort, k_mean, linf_median, k_lwr, k_upr, …)
#> rename: renamed one variable (rate)
#> mutate: new variable 'area' (character) with one unique value and 0% NA
#> filter: removed 363 rows (95%), 19 rows remaining
#> select: dropped 7 variables (cohort, k_mean, linf_median, k_lwr, k_upr, …)
#> rename: renamed one variable (rate)
#> mutate: new variable 'area' (character) with one unique value and 0% NA
#> filter: removed 348 rows (91%), 34 rows remaining
#> select: dropped 7 variables (cohort, k_mean, linf_median, k_lwr, k_upr, …)
#> rename: renamed one variable (rate)
#> mutate: new variable 'area' (character) with one unique value and 0% NA
#> filter: removed 364 rows (95%), 18 rows remaining
#> select: dropped 7 variables (cohort, k_mean, linf_median, k_lwr, k_upr, …)
#> rename: renamed one variable (rate)
#> mutate: new variable 'area' (character) with one unique value and 0% NA
#> filter: removed 367 rows (96%), 15 rows remaining
#> select: dropped 7 variables (cohort, k_mean, linf_median, k_lwr, k_upr, …)
#> rename: renamed one variable (rate)
#> mutate: new variable 'area' (character) with one unique value and 0% NA
#> filter: removed 370 rows (97%), 12 rows remaining
#> select: dropped 7 variables (cohort, k_mean, linf_median, k_lwr, k_upr, …)
#> rename: renamed one variable (rate)
#> mutate: new variable 'area' (character) with one unique value and 0% NA
#> filter: removed 377 rows (99%), 5 rows remaining
#> select: dropped 7 variables (cohort, k_mean, linf_median, k_lwr, k_upr, …)
#> rename: renamed one variable (rate)
#> mutate: new variable 'area' (character) with one unique value and 0% NA

All areas pooled

Code
# One sharpe schoolfield fitted to all data
fit_all <- nls_multstart(
    k_median ~ sharpeschoolhigh_1981(temp = temp, r_tref, e, eh, th, tref = 8),
    data = VBG,
    iter = c(3, 3, 3, 3),
    start_lower = start*0.5,
    start_upper = start*2,
    lower = lower,
    upper = upper,
    supp_errors = 'Y'
    )

# Make predictions on new data
new_data_all <- data.frame(temp = seq(min(VBG$temp),
                                      max(VBG$temp),
                                      length.out = 100))

pred_all <- augment(fit_all, newdata = new_data_all) |> 
  mutate(area = "all")
#> mutate: new variable 'area' (character) with one unique value and 0% NA

# Add t_opt
kb <- 8.62e-05
data.frame(par = names(coef(fit_all)),
           est = coef(fit_all)) |> 
  pivot_wider(names_from = par, values_from = est) |> 
  summarise(t_opt = (eh*th) / (eh + kb*th*log((eh/e)-1)))
#> pivot_wider: reorganized (par, est) into (r_tref, e, eh, th) [was 4x2, now 1x4]
#> summarise: now one row and one column, ungrouped
#> # A tibble: 1 × 1
#>   t_opt
#>   <dbl>
#> 1  19.6

Plot Sharpe Schoolfield fits

Code
# Plot growth coefficients by year and area against mean SST
preds |>
  ggplot(aes(temp, .fitted, color = factor(area, order))) + 
  geom_point(data = VBG, aes(temp, k_median, color = factor(area, order)), size = 0.6) +
  geom_line(aes(temp, .fitted, group = factor(area)), linewidth = 1) +
  geom_line(data = pred_all, aes(temp, .fitted),
            linewidth = 1, inherit.aes = FALSE, linetype = 2) +
  #scale_color_manual(values = colors, name = "Area") +
  scale_color_viridis(name = "Area", direction = -1, discrete = TRUE) +
  theme(plot.title = element_text(size = 15, face = "bold")) +
  theme(axis.text.x = element_text(angle = 0)) +
  theme(axis.text = element_text(size = 12), axis.title = element_text(size = 15)) +
  guides(color = guide_legend(override.aes = list(size = 1))) +
  scale_x_continuous(breaks = seq(-5, 20, 1)) +
  labs(x = "Temperature (C)", y = "Median von Bertalanffy growth coefficient, k") +
  NULL

Code

ggsave(paste0(home, "/figures/sharpe_school.pdf" ), width = 17, height = 17, units = "cm")

Test plotting Linf aswell…

By area

Code
model <- 'sharpeschoolhigh_1981'

# Get starting values on full dataset for Sharpe-Schoolfield model
dat <- VBG |>
  select(linf_median, temp) |>
  rename(rate = linf_median)
#> select: dropped 8 variables (cohort, area, k_mean, k_median, k_lwr, …)
#> rename: renamed one variable (rate)

lower <- get_lower_lims(dat$temp, dat$rate, model_name = model)
upper <- get_upper_lims(dat$temp, dat$rate, model_name = model)
start <- get_start_vals(dat$temp, dat$rate, model_name = model)
  
# Sharpe-Schoolfield model fit to data for each area
preds <- NULL

for(a in unique(VBG$area)) {
  
  # Get data
  dat <- VBG |> 
    filter(area == a) |> 
    select(linf_median, temp, area) |> 
    rename(rate = linf_median)
  
  # Fit model
  fit <- nls_multstart(
    rate ~ sharpeschoolhigh_1981(temp = temp, r_tref, e, eh, th, tref = 8),
    data = dat,
    iter = c(3, 3, 3, 3),
    start_lower = start*0.5,
    start_upper = start*2,
    lower = lower,
    upper = upper,
    supp_errors = 'Y'
    )
  
  # Make predictions on new data
  new_data <- data.frame(temp = seq(min(dat$temp), max(dat$temp), length.out = 100))
  
  pred <- augment(fit, newdata = new_data) |>
    mutate(area = a)
  
  # Add to general data frame
  preds <- data.frame(rbind(preds, pred))
  
}
#> filter: removed 319 rows (84%), 63 rows remaining
#> select: dropped 7 variables (cohort, k_mean, k_median, k_lwr, k_upr, …)
#> rename: renamed one variable (rate)
#> mutate: new variable 'area' (character) with one unique value and 0% NA
#> filter: removed 328 rows (86%), 54 rows remaining
#> select: dropped 7 variables (cohort, k_mean, k_median, k_lwr, k_upr, …)
#> rename: renamed one variable (rate)
#> mutate: new variable 'area' (character) with one unique value and 0% NA
#> filter: removed 341 rows (89%), 41 rows remaining
#> select: dropped 7 variables (cohort, k_mean, k_median, k_lwr, k_upr, …)
#> rename: renamed one variable (rate)
#> mutate: new variable 'area' (character) with one unique value and 0% NA
#> filter: removed 347 rows (91%), 35 rows remaining
#> select: dropped 7 variables (cohort, k_mean, k_median, k_lwr, k_upr, …)
#> rename: renamed one variable (rate)
#> mutate: new variable 'area' (character) with one unique value and 0% NA
#> filter: removed 334 rows (87%), 48 rows remaining
#> select: dropped 7 variables (cohort, k_mean, k_median, k_lwr, k_upr, …)
#> rename: renamed one variable (rate)
#> mutate: new variable 'area' (character) with one unique value and 0% NA
#> filter: removed 344 rows (90%), 38 rows remaining
#> select: dropped 7 variables (cohort, k_mean, k_median, k_lwr, k_upr, …)
#> rename: renamed one variable (rate)
#> mutate: new variable 'area' (character) with one unique value and 0% NA
#> filter: removed 363 rows (95%), 19 rows remaining
#> select: dropped 7 variables (cohort, k_mean, k_median, k_lwr, k_upr, …)
#> rename: renamed one variable (rate)
#> mutate: new variable 'area' (character) with one unique value and 0% NA
#> filter: removed 348 rows (91%), 34 rows remaining
#> select: dropped 7 variables (cohort, k_mean, k_median, k_lwr, k_upr, …)
#> rename: renamed one variable (rate)
#> mutate: new variable 'area' (character) with one unique value and 0% NA
#> filter: removed 364 rows (95%), 18 rows remaining
#> select: dropped 7 variables (cohort, k_mean, k_median, k_lwr, k_upr, …)
#> rename: renamed one variable (rate)
#> mutate: new variable 'area' (character) with one unique value and 0% NA
#> filter: removed 367 rows (96%), 15 rows remaining
#> select: dropped 7 variables (cohort, k_mean, k_median, k_lwr, k_upr, …)
#> rename: renamed one variable (rate)
#> mutate: new variable 'area' (character) with one unique value and 0% NA
#> filter: removed 370 rows (97%), 12 rows remaining
#> select: dropped 7 variables (cohort, k_mean, k_median, k_lwr, k_upr, …)
#> rename: renamed one variable (rate)
#> mutate: new variable 'area' (character) with one unique value and 0% NA
#> filter: removed 377 rows (99%), 5 rows remaining
#> select: dropped 7 variables (cohort, k_mean, k_median, k_lwr, k_upr, …)
#> rename: renamed one variable (rate)
#> mutate: new variable 'area' (character) with one unique value and 0% NA

All areas pooled

Code
# One sharpe schoolfield fitted to all data
fit_all <- nls_multstart(
    linf_median ~ sharpeschoolhigh_1981(temp = temp, r_tref, e, eh, th, tref = 8),
    data = VBG,
    iter = c(3, 3, 3, 3),
    start_lower = start*0.5,
    start_upper = start*2,
    lower = lower,
    upper = upper,
    supp_errors = 'Y'
    )

# Make predictions on new data
new_data_all <- data.frame(temp = seq(min(VBG$temp),
                                      max(VBG$temp),
                                      length.out = 100))

pred_all <- augment(fit_all, newdata = new_data_all) |> 
  mutate(area = "all")
#> mutate: new variable 'area' (character) with one unique value and 0% NA

# Add t_opt
# kb <- 8.62e-05
# data.frame(par = names(coef(fit_all)),
#            est = coef(fit_all)) |> 
#   pivot_wider(names_from = par, values_from = est) |> 
#   summarise(t_opt = (eh*th) / (eh + kb*th*log((eh/e)-1)))

Plot Sharpe Schoolfield fits

Code
# Nice color palette, order based on average temperature!
# Plot growth coefficients by year and area against mean SST
preds |>
  ggplot(aes(temp, .fitted, color = factor(area, order))) + 
  geom_point(data = VBG, aes(temp, linf_median, color = factor(area, order)), size = 0.6) +
  geom_line(aes(temp, .fitted, group = factor(area)), linewidth = 1) +
  geom_line(data = pred_all, aes(temp, .fitted),
            linewidth = 1, inherit.aes = FALSE, linetype = 2) +
  # geom_smooth(aes(temp, .fitted, color = factor(area, order)), 
  #             method = "gam", formula = y ~ s(x, k = 4)) +
  # geom_smooth(aes(temp, .fitted), inherit.aes = FALSE,
  #             method = "gam", formula = y ~ s(x, k = 4), linetype = 2) +
  #scale_color_manual(values = colors, name = "Area") +
  scale_color_viridis(name = "Area", direction = -1, discrete = TRUE) +
  theme(plot.title = element_text(size = 15, face = "bold")) +
  theme(axis.text.x = element_text(angle = 0)) +
  theme(axis.text = element_text(size = 12), axis.title = element_text(size = 15)) +
  guides(color = guide_legend(override.aes = list(size = 1))) +
  scale_x_continuous(breaks = seq(-5, 20, 1)) +
  labs(x = "Temperature (C)", y = "Asymptotic size (mm)") +
  NULL

Code
ggplot(VBG, aes(linf_median, k_median)) + 
  geom_point()